Deposition of magnetic particles: A computer simulation study 
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We report a Monte Carlo simulation of deposition of magnetic particles on a one-dimensional 
substrate. Incoming particles interact with those that are already part of the deposit via a dipole- 
dipole potential. The strength of the dipolar interaction is controlled by an effective temperature 
T*, the case of pure diffusion- limited deposition being recovered in the limit T* — » oo. Preliminary 
results suggest that the fractal dimension of the deposits does not change with temperature but that 
there is a (temperature-dependent) cross-over from regimes of temperature-dependent to universal 
behaviour. Furthermore, it was found that dipoles tend to align with the local direction of growth. 
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I. INTRODUCTION 



Magnetic particles are a key ingredient of many mod- 
ern data recording and storage devices, from music 
tapes to computer hard disks Q. For these applica- 
tions smooth, regular magnetic layers are usually desired. 
However, dipoles, both magnetic and electric, display a 
fondness for arranging themselves into highly inhomo- 
geneous structures. This is a consequence of the very 
strong anisotropy of the dipole-dipole interaction, which 
couples the orientations of the dipole moments with that 
of the interparticle vector. Because the potential energy 
at a fixed separation is lowest for a head-to-tail geome- 
try, chain formation is particularly favoured in ferrofluids 
(dispersions of ferromagnetic particles) || or electrorhe- 
ological fluids (dipersions of highly polarisable particles 
in solvents with low dielectric constant) || in magnetic 
and electric fields, respectively. Whether such chaining 
can occur in zero field in the absence of any interactions 
other than dipolar, is experimentally uncertain; it has 
been seen in simulations, but is not yet fully accounted 
for theoretically. Likewise, what is the true equilibrium 
structure of a solid of hard magnetic particles is still un- 
settled (see, e.g., j| and references therein). 

It is therefore of great importance, from a practical 
as well as from a fundamental point of view, to inves- 
tigate the influence of true long-range magnetic dipolar 
forces on the geometry of particle aggregates, so as to be 
able to exert better control over them. This is especially 
relevant to the very novel field of self-assembled nanos- 
tructured magnetic materials, where the aim is to allow 



different microscopic components to organise themselves 
into complex functional patterns once their interactions 
have been appropriately tailored Many of these 

devices, either existing or at the design stage, have low 
dimensionalities (e.g., wires and films), at which simula- 
tions of model systems have revealed the chaining ten- 
dency to be particularly strong |||7]|y|. 

To our knowledge, there is detailed only one study of 
how dipolar interactions alone affect growth. Pastor- 
Satorras and Rubf simulated an off-lattice, two- 
dimensional particle-cluster aggregation model. They 
found a monotonic variation of the fractal dimension 
of the aggregates as a function of temperature (i.e., 
dipole strength), the limit of diffusion- limited aggrega- 
tion (DLA) being recovered at high temperatures. In ad- 
dition, a separate investigation by the same authors has 
shown that highly structured layers could be obtained 
at low temperatures by sequential adsorption of dipolar 
particles . See |l2|,|l3| for related work on other sys- 
tems. 

Here we report on a simulation of dipole deposition 
on a one-dimensional (Id) substrate (i.e., a line). In the 
limit of zero magnetic moment this reduces to diffusion- 
limited deposition (DLD), which should exhibit the same 
geometrical properties as DLA JlJ] ; our work thus serves 
as a both a check and an extension of Pastor-Satorras 
and Rubf's to the case of an infinite (in one spatial di- 
mension) system. In particular, we want to ascertain: 
firstly, whether the fractal dimension actually changes 
owing to the strong anisotropy of the dipolar interaction: 
and, secondly, what is the correlation between the orien- 
tations of dipoles in the aggregate and its direction(s) of 
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growth. For computational convenience our dipoles are 
restricted to reside at the sites of a two-dimensional (2d) 
square lattice (although they can point in any direction 
of three-dimensional (3d) space): we assume that any ef- 
fects coming from this discretisation of space are much 
smaller than those of the interparticle potential. That 
this should be so is suggested by results for DLA 
(but remains of course to be confirmed by full off-lattice 
calculations). Furthermore, our analysis in terms of the 
concepts of fractal geometry presumes that all our de- 
posits are self-similar over some lengthscale larger than 
the mesh spacing but smaller than the deposit size jL4| ; 
again, this need not be true of the smallest deposits, but 
these contain only a very small fraction of the total num- 
ber of particles. 

The present paper is a natural progression to non- 
equilibrium processes from our earlier researches on the 
thermodynamics and phase equilibria of dipolar fluids [Q . 
It is organised as follows: in section || we describe our 
model and the simulation method employed. Then in 
section [II we present and discuss our results, specifically 



comparing them with those of Pastor-Satorras and Rubf 
[0. Section |v| summarises our findings and outlines 
prospects for future research. Technical details pertain- 
ing to the treatment of the long range of the dipole-dipole 
interaction are relegated to two appendices. 



II. MODEL AND SIMULATION METHOD 

Our simulations were performed on a (1+1)- 
dimensional square lattice of width L — 800a sites and 
any height that can accommodate N dipoles, where a 
is the mesh spacing and the adsorbing substrate coin- 
cides with the bottom row (henceforth we take a = 1). 
Periodic boundary conditions are imposed in the direc- 
tion parallel to the substrate. Each particle carries a 3d 
dipole moment of strength [i and interact through the 
pair potential 



>D 



(1) 2 ) = g - • ha)^ • ?i 2 ) --Mi > (!) 



' 12 



where ri2(> a) is the distance between particles 1 and 
2, f 12 is the two-dimensional (2d) unit vector along the 
interparticle axis, and -p^ and p^ are the 3d unit vec- 
tors in the direction of the dipole moments of particles 1 
and 2 respectively. Finally, '1' and '2' denote the full set 
of positional and orientational coordinates of particles 1 
and 2. 

A particle is introduced at a lattice site (xi n , H max + 
AL), where Xi n is a random integer in the interval [1, L], 
H m ax is the maximum distance from the substrate to 
any particle in the deposit, and A is a constant. The 
dipole moment of the released particle is oriented at ran- 
dom. The particle then undergoes a random walk by a 
series of jumps to nearest-neighbour lattice sites, while 
experiencing dipolar interactions with the particles that 



are already attached to the deposit. We incorporate the 
effects of these interactions through a Metropolis algo- 
rithm. If the deposit contains N particles, then the in- 
teraction energy of the (N + l)th incoming particle (the 
random walker) with the particles in the deposit is given 

by E(r,-p) = EiIiEn^(*' iV + l )> where r and ^ 
are the current position and the dipole orientation of the 
random walker respectively (r is a 2d vector). Then we 
randomly choose a new position r' (|r — r'| = a) and 
a new dipole orientation pi for the random walker; this 
displacement is accepted with probability 



p = min < 1 , exp 



£(rV) - E(r,*,) 



(2) 



T* = ksTa 3 /^ 2 is an effective temperature, inversely 
proportional to the dipolar energy scale. In the limit 
T* — > only displacements that lower the energy 
E(r' ,-pJ) are accepted. On the other hand, in the limit 
T* — ► oo all displacements are accepted and our model 
reduces to the well-known DLD jl4|] . 

The long range of the dipole-dipole interaction was 
treated by the Ewald sum method (see Appendix A for 
details). In our simulations we set a = 10/L, for which it 
suffices to re tain terms with n = in the real space sum 
of equation (All). The sum in reciprocal space extends 
over all lattice points k = 2irn/L with |n| < 16, whereas 
the sum in real space is truncated at L/2. 

The particle eventually either contacts the deposit (i.e., 
becomes a nearest neighbour of another particle that is 
already part of the deposit) or moves away from the sub- 
strate. In the latter case, if the particle reaches a dis- 
tance from the substrate greater than H max + 2AL, it is 
removed and a new one is launched. Once a particle has 
reached the substrate or the deposit, its dipole relaxes 
along the direction of the local field created by all other 
particles in the deposit. In all simulations reported here 
we took A = 1; larger values of A were tested and found 
to give the same results, but with drastically increased 
computation times. 



III. RESULTS AND DISCUSSION 



Four effective temperatures were considered: T* = 
10" 1 (28 deposits), 10" 2 (41 deposits), 10" 3 (42 deposits) 
and 10~ 4 (54 deposits), chosen to be in the range where 
the fractal dimension of dipolar DLA clusters is expected 
to change jHJ. Each deposit contains 50000 dipoles. We 
have also generated 30 DLD deposits on the same lattice 
by this same method, with T* — oo; known results for 
DLD (see, e.g., p6[|) were used to check the validity of 
our algorithm. On the other hand, comparison between 
these and the results for finite temperatures reveal the 
effect of dipolar interactions on DLD. 
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P (h) = ( I Y / v(hh)), 
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FIG. 1. Snapshots of two deposits obtained for T* = 10 1 
(black) and T* = 10" 4 (grey). 



Figure |l| is a snaphsot of two deposits for T* = 10 _1 
(black) and T* = 10~ 4 (grey, only part of the deposit is 
shown - in fact, it grows up to a height of about 8000). 
Both deposits have the same general appearance, already 
observed in DLD: they consist of several trees compet- 
ing to grow. As the size of the deposit increases, fewer 
and fewer trees 'survive' (i.e., carry on growing), as a 
consequence of the so-called shielding or screening effect. 
From figure this seems more pronounced at lower tem- 
peratures, since above a height of 1000 (about 1/8 of the 
maximum height attained by this deposit) only one tree 
survives. 

In order to compare quantitatively the results obtained 
for different temperatures we have measured the mean 
density of dipoles p{h) at a height h, 





FIG. 2. (a) Mean density p(h) of deposits at height h: ran- 
dom deposition (solid line), T* = 10 _1 (dotted line), 10~ 2 
(dashed line), 10~ 3 (long-dashed line) and 10 -4 (dot-dashed 
line). (b) h* p(h) vs h/h*: random deposition (circles), 
T* = 1CF 1 (squares), 10~ 2 (diamonds), 10~ 3 (triangles up) 
and 1CF 4 (triangles down). The solid line is a linear fit to the 
DLD results in the range of the graph; it has slope —0.28. 

where r)(i,h) is 1 (0) if the site with coordinates (i, h) 
is occupied (unoccupied), and the average (denoted by 
angular brackets) is taken over all available deposits at 
each temperature. p(h) is plotted in figure 0a: all the 
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curves have similar shapes, with a smooth decrease at 
small h, levelling off (saturating) at intermediate h, and 
with a sharp drop at large h, when the top of the de- 
posit is reached. It is immediately noticed that the finite- 
temperature curves differ from that for DLD in one im- 
portant respect: the density at a given height and the 
maximum height h* attained by the deposits vary with 
temperature. These variations are monotonic: p(h) is 
smaller and h* is larger at lower temperatures and thus 
the increase in the strength of dipolar interactions en- 
hances the shielding effect. 

For DLD, p{h) was found to be of the form E3|, 



p{h) oc h- a g(h/h*), 



(4) 



where a, the so-called codimensionality, is the difference 
between the dimension of space d (= 2, in the presente 
case) and the fractal dimension D of the deposit, h* is 
the maximum height, and g{x) is a universal function. 
g{x) w 1 when a; < 1 and decays faster than any power 
of x when x — > 1. In order to compare DLD and finite- 
temperature results, we propose a general form for p(h), 
inspired by equation (Q): 



p(KT*) = A(T*)h~ a g(h/h*), 



(5) 



where A(T*) is some unknown function of T* only. It is 
easily seen that A(T*) can be found as a function of h* 
and of the number of particles in the deposit N, by using 
the normalization condition. 



h* 

v = l x Vp(A,r) 



h=l 



L 



p{h,T*)dh, 



whence 



A(T*) 



L fi/h* x a 9(x)dx 



(6) 



(7) 



Since g(x) is a universal function and L and N are the 
same for all the deposits we are analysing, equation (|^) 
becomes 



p(h, T*)h* cx 



g(h/h*). 



(8) 



In figure ||b we plot h*p(h,T*) as a function of h/h*, 
for h/h* <C 1 on a log-log scale. The data points for 
DLD follow, as expected, a straight line. A linear re- 
gression gives a w 0.27 and thus D w 1.73, in good 
agreement with what was obtained previously by several 
other methods |l7]Jltj ] . If the full functional dependence 
of p{h,T*) were captured by equation (||), two possibil- 
ities would arise: (i) a is T*-independent and all curves 
are paralell straight lines; (ii) a depends on T and all 
curves are straight lines with different slopes. However, 
we arrive at neither of these scenarios, so the situation 
is a little more complex. If we were to interpret our re- 
sults in terms of a function similar to equation (H), then 



g(x) would also need to have an explicit temperature de- 
pendence. This dependence should be able to describe 
the trends observed in figure ||b for finite temperatures: 
a crossover between an approximately linear regime for 
very low relative heights, characterised by an exponent 
greater than a = 0.27; and a linear regime for intermedi- 
ate heights, characterised by roughly the same exponent 
as DLD. 

In jlTJ] it was argued that the change in the confor- 
mational properties of DLA clusters introduced by the 
presence of dipolar interactions could be interpreted as 
a change in their fractal dimension. The fractal dimen- 
sion for each temperature was determined by measur- 
ing the dependence of the radius of gyration of dipolar 
DLA clusters on the number of particles in a cluster. 
Between T* = 10" 1 and T* = 10~ 4 a fractal dimen- 
sion was obtained ranging from about 1.7 to about 1.2. 
We have performed linear fits to the data shown in fig- 
ure |^b, using only those points corresponding to heights 
below the region where the crossover referred to above 
seems to take place. These points follow straight lines 
with temperature-dependent slopes ranging from 0.3 (for 
T* = 1Q- 1 ) to 0.6 (for T* = 10~ 4 ). On the basis of 
the analysis of just this part of the deposit, we obtain 
an (apparent) variation of the fractal dimension of the 
deposits with T*, from D = 1.7 to D = 1.4. We con- 
clude, as was already pointed out in and seems to 
be confirmed by the present work, that the results of 
|]l0f can be interpreted in terms of a crossover between a 
temperature-dependent fractal dimension at short length 
scales, tofl« 1.7 at long length scales, with a crossover 
height that itself depends on temperature. This conclu- 
sion must, however, be tested with longer simulations at 
the lowest temperatures, where the statistics are a little 
poorer. 

There are other routes to estimating the fractal di- 
mension of the deposits. We shall use one other to show 
that it is not necessary to assume that the finite tem- 
perature deposits have a fractal dimension different from 
DLD. The mean height of the upper surface, h m . when 
the deposit contains M particles, is defined as |lq| , 



h m (M) = {jJ2h max (i,M)), 



(9) 



where h max (i,M) is the maximum height of the occu- 
pied sites of column i when there are M particles in the 
deposit. In a DLD deposit this quantity is expected to 
scale with M, as flq] 



(10) 
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FIG. 3. (a) Mean height of the upper surface, h m , as a 
function of the number of particles, M. The lines are as in 
figure ^a. (b) Blowup of the large-M region. 

The exponent <f> is related to the codimcnsionality 
a = d — D by (/) — j^— and to the fractal dimension 
by D = d — 1 + In figure [§] we plot our results for 
h m {M). The scaling law, equation (|9f), is known to be 
valid in the limit M — > 00 [p.6| . As in ||16| , we have per- 
formed several linear regressions for large M, in the range 
Mi < M < M 2 , for (Mx,M 2 ) = (0.57V, N), (0.25N,N), 
(0.1N,N) and (0.257V, 0.5 JV) (recall that N = 50000). 
For every temperature and every range considered we 



found that 1.33 < 4> < 1.44, which corresponds to a frac- 
tal dimension, 1.69 < D < 1.75. Moreover, we have 
found no evidence of any regular variation of <j> with T* 
over a given range of M. Figure ||a shows that, in the 
initial stages of growth, the mean height of the upper 
surface grows identically at every temperature. There is 
then a crossover region at intermediate stages when the 
less dense deposits grow slightly faster. Finally, in the 
later stages all the deposits grow at the same rate regard- 
less of temperature, as evidenced by figure ||b. However, 
note that, as is clear from figure ||a, if the deposits had 
been allowed to grow only to intermediate stages (e.g., up 
to M = 10000), an increasing value of <fi with increasing 
temperature would have obtained, and thus an apparent 
variation of D with T* , with the same trends as observed 
through the calculation of p(h). 

We conclude this analysis by attempting to make a 
first connection between the orientation of the dipoles in 
the deposit and its growth. Figure |J is a snapshot of part 
of a deposit for T* = 10 _1 . Dipoles whose horizontal (or 
lateral) component is smaller (greater) in absolute value 
than their vertical component are coloured black (grey). 
Since we have verified that the z-component (i.e., out-of- 
plane) of the dipoles in the deposit is zero after a short 
time, figure [| suggests that the dipoles tend to align with 
the direction of growth of the deposit at the site where 




FIG. 4. Detail of a deposit for T* = 1CT 1 . Sites whose 
dipoles make an angle of absolute value smaller (larger) than 
7r/4 with the vertical axis are shown in black (grey). 

In order to make this idea more quantitative, we have 
measured the angles lo between the direction of the dipole 
moments of all incoming particles and the direction of 
growth at their point of attachment to the deposit. We 
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have done so by recording whether a new dipole becomes 
attached to the substrate due to a neighbour positioned 
to its left or to its right (lateral growth: the relevant an- 
gle is that between the dipole and the horizontal axis), 
or above or below it (vertical growth: the relevant an- 
gle is now that between the dipole and the vertical axis) . 
We did not take into account particles that attach to the 
deposit having both vertical and horizontal neighbours. 
Once these angles were collected, for every deposit at 
each temperature, we constructed a frequency histogram 
by dividing the interval (— tt,tt) into 1000 sub-intervals. 
In figure |5| we plot these results for lateral and vertical 
growth at T* = 10" 1 and T* = 10~ 4 . 




b) 



| 0.02 
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FIG. 5. Frequency histogram of the angle to between 
dipole orientation and the local direction of growth, for (a) 
T* = 10" 1 and (b) T* = 10" 4 . The solid lines correspond to 
lateral growth and the dotted lines to vertical growth. 

Because all curves have period 7r and are even, only 
the interval (— 7r/2,7r/2) is shown. The main feature of 
all the curves are the strong peaks at u = 0, imply- 
ing that most dipoles align in the direction of growth of 



the deposit. This is a consequence of the fact that the 
lowest-energy configuration of two dipoles a fixed dis- 
tance apart is head-to-tail along the direction of the in- 
terdipole vector. These peaks are more pronounced at 
the lower temperature, and the peak for vertical growth 
is a little higher than that for lateral growth at both tem- 
peratures, implying that growth in the vertical direction 
is more likely to happen with dipoles aligned vertically 
than growth in the lateral direction with dipoles aligned 
horizontally. 

The curves exhibit some other, lower, peaks. There 
is a broad, low peak around lu = ir/2 at T* = 10 _1 , 
which corresponds to lateral growth with vertically- 
aligned dipoles ('black' horizontal branches in figure ||), 
or to vertical growth with horizontally-aligned dipoles 
('grey' vertical branches in figure ^), at high tempera- 
tures. This can be explained by noting that the second- 
lowest minimum of the interaction energy at fixed separa- 
tion is for two antiparallel dipoles. These peaks seem to 
disappear at T* — 10~ 4 , suggesting that energetic effects 
become more important as the temperature is lowered. 

There are several other peaks, occurring at the same 
angles for both lateral and vertical growth, whose posi- 
tions seem unaffected by changing the temperature. At 
the present stage of our work, we can only speculate as 
to their origin. We believe these peaks come from a com- 
bination of lattice effects and the properties of the dipo- 
lar interaction. It is actually known that the minimum- 
energy arrangement of n(> 3) dipoles is obtained by plac- 
ing them at the vertices of a regular n-sided polygon, 
tangent to the circumscribing circle. Thus four dipoles 
on a square lattice will minimse their energy by making 
7r/4 angles with the horizontal and vertical axes, which 
might explain the peaks observed in figure || at that an- 
gle. The remaining peaks may likewise correspond to 
other arrangements of dipoles realising other minima of 
the energy of sets of dipoles on a square lattice. 



IV. CONCLUDING REMARKS 

We have simulated the deposition of dipolar particles 
on a Id substrate using a lattice model. Our findings sug- 
gest that the fractal dimension of the deposits is the same 
as for DLD and hence unaffected by the dipolar interac- 
tions, but also that there is a crossover from temperature- 
dependent to temperature independent behaviour which 
can be very broad. A fuller characterisation in terms of 
the height-height correlation function, tree size distribu- 
tions and density scaling with system size, is in progress 
and will be published elsewhere. These quantities would 
provide additional routes to the fractal dimension, thus 
allowing us further to verify (or falsify) our preliminary 
conclusions. Growth and roughness exponents will also 
be calculated. 

We have now started work on the off-lattice version 
of the present model, so as to be free from any possible 
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artifacts arising from the discretisation of space. Results 
so far suggest that the lattice has a very small effect, but 
we are currently somewhat limited by the very high com- 
putational cost of evaluating the interactions between a 
particle and all its periodic images at every step. More 
efficient algorithms are being developed to tackle this. 
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APPENDIX A: EWALD SUMS FOR THE 
DEPOSITION PROCESS 



We have generalised to arbitrary dimensions a method 
proposed by Grzybowski et al. |l8j for evaluating Ewald 
sums. Here for simplicity we just present results for our 
case d = 1 (where d is the number of dimensions in which 
we impose periodic boundary conditions). Our simula- 
tion box consists of a rectangle with a base of length L 
and any height that can accommodate N dipoles. The 
dipoles are always three-dimensional vectors, with po- 
sition vectors R^, i = 1, 2, . . . , TV. The simulation box is 
repeated in the x (horizontal) direction, giving rise to a 
regular lattice whose sites are located at n = (n, 0)L. Let 
fj, and R be the dipole moment and position of an incom- 
ing particle. Then, the distance between the incoming 
particle in the origin cell and another in an image cell is 
Ti = R (Rj + n), i = 1 , 2, . . . , N. The total interaction 
energy between the incoming particle and the N particles 
in the box and their infinite replicas is 



\H- (r. ( +n)][n i ■ (r;+n)] 



(Al) 

According to the geometry of the system = [xi , t/, , 0) 
where Xi (yi) is the horizontal (vertical) distance between 
the incoming particle and a particle in the deposit. Note 
that the incoming particle does not interact with its own 
images. 

Introducing the notation 




' n + rl d 



r, c 



e -ic-(n+r) 
Z-^ n 4- r 5 ' ' 



(A2) 
(A3) 



allows us to express the total energy as 

N 
i=l 

N 

+ 3^ ( M • V c ) (m< ' V c )0(r*,c)| c=o . (A4) 

i=l 

To calculate ip{r) and #(r, c) we use the identities: 

1 



x~ 2u = 



r(«) Jo 



t u - L e- xt dt, (A5) 



e -t|r+n| 2 -ic-(r+n) 



7T 
tl? 



1/2 



£k-r 



(k + c) 2 \ 



k = 27r(A,0)/L with k integer is a reciprocal-lattice vec- 
tor and r and n are as above. Equation flAg ) is a direct 
consequence of the definition of the Gamma function, 
while equation ( |A6| ) is a form of the Poisson summation 
formula for d = 1, which is the dimensionality of the 
lattice formed by repeating the box. 
For ip(r) we set u = 3/2, leading to 



m = 



(A7) 



The sum over direct-lattice vectors converges fast for 
large t, while that over reciprocal-lattice vectors does so 
for small t. We therefore choose an arbitrary separation 
parameter a 2 for the t integration and decompose the 
lattice sum into two terms: 

a 2 

+ t l/2 e -t\r+^ dL 



(A8) 



Taking into account that — i|r+n| 2 = — t\x+n\ 2 ~ ty 2 and 
using the Poisson summation formula, equation ( |A6| ) , we 
arrive at 

9 /*°° 

In the case of 9(r, c) we set u = 5/2 with the result 
"(r^) = w^E e ^' (r+n) J a2 t 3/2 e- tlr+ " l2 dt (AfO) 



r(5/2) 
4 

+ 3LE 



ikx—iycy 



texp( -ty z + —\dt. 
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We now substitute ?/>(r) and 0(r, c) into equation (A4) 
for E: 

£ 4^{ (,i '* (a, ' J) 

* i n 

- 2 ^ • (Ti + n) /x, • (r 4 + n) 7 3/2 (a, /3)| 



L 



i k 



7 Mj/^^ etfcx<j i( a ' fc ) 

t k 

■r^^MsMixfeV^J-^a.j/i,*;) (All) 

i k^O 

y^y^ k(u x ix iy + ti lx ^y^y i e lkx 'Jo(a 1 y l ,k), 



i k/0 



where /? = |r^ + n| and the integrals are given by (see 
Appendix B for details) 



h 



/2 



'3/2 



(a,\ri +n 



(a, \ri + n 



J v (a,yi,k) = 



Wi + n • 



^ I erfc(a|r i : + „|), (A12) 



|r 4 + 
3v^ 



4|r, + n| 5 

t v e -ty1-k*IU di 



2|r i + n| 2 
erfc( a\ti + n 



(A13) 



Since E is real, equation ( All ) can be further simpli- 
fied by replacing exp(ikxi) by cos(focj) + i sin(fcEi). Note 
also that it is even in k and that the case k = can be 
treated analytically. Therefore, no distinction is made 
in the code between simulation box images with positive 
and negative fc, and the case k = is considered sepa- 
rately. 



APPENDIX B: EVALUATION OF SOME 
RELEVANT INTEGRALS 

The following definitions and results will be useful: 

(Bl) 
(B2) 



e~ vt dt = -«/-, 



1 TT 

2V^ 



erf(w) = —= I e dt, 
V 71 " Jo 

2 

erfc(u) = 1 — erf (it) = —p 



(B4) 
(B5) 



We are interested in 7l « and 7 3 m . They can be evaluated 
with the help of 7(— 1/2), which is easy: 



I„(a,/3) = I t v e~ tl3 dt, 
J v (a,a,b) = I t v e~ at - b/t dt. 



/-i/ 2 («,/3) = ^erfc(V) 



(B6) 



where we have made the change of variable t — z 2 //3 2 . 
Now the cases v — 1/2 and v = 3/2 can be easily found 
by integrating by parts: 



h/2(a,P) 
73/2 («>/3) 



ae 



+ 



2 2fr 



5 erfc M> 



(B7) 



-0 2 a 2 



ft 2 



30? 



a + W) + W e " ic[af3)l (B8) 



having resorted to the changes of variables u = t 1 / 2 and 
u' = —e~P *. Turning next to J v , the values of v we are 
interested in depend on the system dimensionality. For 
d = 1, v = —1,0,1 and no analytical results are avail- 
able, in contrast with the d — 2 case. One of the three 
integrals, however, can be expressed in terms of the other 
two: 

Y J-i + Jo - y 2 Ji = a 2 cxp ( - y 2 o? - ^) . (B9) 



e _t dt. (B3) 



Two classes of integrals need to be performed: 
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